##### functions to be used in deterministic satellite-debris model dynamic programming solver script

# Collision rate 
L <- function(S,D,params) {
	aSS <- params$aSS
	aSD <- params$aSD
	1-exp( -(aSS*S) -(aSD*D))
}

# Derivatives of collision rate
L_S <- function(S,D,params) {
	aSS <- params$aSS
	aSD <- params$aSD
	aSS*exp(- aSS*S - aSD*D)
}

L_D <- function(S,D,params) {
	aSS <- params$aSS
	aSD <- params$aSD	
	aSD*exp(- aSS*S - aSD*D)
}

# debris growth function
G <- function(S,D,params) {
	aSS <- params$aSS
	aSD <- params$aSD
	aDD <- params$aDD
	bSS <- params$bSS
	bSD <- params$bSD
	bDD <- params$bDD
	bSS*(1-exp(-aSS*S))*S + bSD*(1-exp(-aSD*D))*S + bDD*(1-exp(-aDD*D))*D
}

# average new fragments function S derivative
G_S <- function(S,D,params) {
	aSS <- params$aSS
	aSD <- params$aSD
	bSS <- params$bSS
	bSD <- params$bSD
	bSS*(1-exp(-aSS*S)) + bSS*aSS*S*exp(-aSS*S) + bSD*(1-exp(-aSD*D))
}

# average new fragments function D derivative
G_D <- function(S,D,params) {
	aSD <- params$aSD
	aDD <- params$aDD
	bSD <- params$bSD
	bDD <- params$bDD
	aSD*bSD*S*exp(-aSD*D) + bDD*(1-exp(-aDD*D)) + aDD*bDD*exp(-aDD*D)*D
}

# Satellite law of motion 
S_ <- function(X,S,D,params) {
	S*(1 - L(S,D,params)) + X
}

# Debris law of motion
D_ <- function(X,S,D,params) {
	stock <- D*(1-d) + G(S,D,params) + m*X
	stock[which(stock=="NaN")] <- D
	stock[which(stock=="NA")] <- D
	stock
}

# open access equilibrium condition
eqmcond <- function(X,S,D,params) {
	fe_eqm <- params$fe_eqm
	S_next <- S_(X,S,D,params) 
	D_next <- D_(X,S,D,params)
	L(S_next,D_next,params) - fe_eqm
}

# analytical expression for launch rate assuming negative exponential collision rate
eqm_launch_rate <- function(S,D,params) {
	fe_eqm <- params$fe_eqm
	num <- log(1-fe_eqm) + aSS*(1-L(S,D,params))*S + aSD*(D*(1-d)+G(S,D,params))
	den <- aSS + aSD*m
	launches <- -num/den
	launches <- pmax(launches,0)
	return(launches)
}

# analytical expression for open-access steady-state launch rate assuming negative exponential collision rate
oass_launch_rate <- function(S,D,params) {
	fe_eqm <- params$fe_eqm
	num <- log(1-fe_eqm) + aSS*S + aSD*D
	den <- aSS + aSD*m
	launches <- -num/den
	launches <- pmax(launches,0)
	return(launches)
}

# analytical expression for the equilibrium steady-state satellite stock as a function of the debris stock assuming negative exponential collision rate
eqm_sat_stock <- function(D,params) {
	fe_eqm <- params$fe_eqm
	value <- -(1/aSS)*( log(1-fe_eqm) + aSD*D )
	value <- pmax(0,value)
	return(value)
}

# One-period fleet returns
one_p_return <- function(X,S,params) {
	p <- params$p
	F <- params$F
	p*S - F*X
}

# Open-access steady state condition with constant launch rate
ss_cond_const <- function(init,params,constant_launch_rate) {
	S <- init[1]
	D <- init[2]
	X <- constant_launch_rate
	S_imbalance <- S_(X,S,D,params) - S
	D_imbalance <- D_(X,S,D,params) - D
	return(c(S_imbalance,D_imbalance))
}

# Open-access steady state condition with general launch rate
ss_cond <- function(init,params) {
	S <- init[1]
	D <- init[2]
	X <- eqm_launch_rate(S,D,params)
	S_imbalance <- S_(X,S,D,params) - S
	D_imbalance <- D_(X,S,D,params) - D
	return(c(S_imbalance,D_imbalance))
}

# Optimal steady state condition
ss_cond_optimal <-function(init,params){
	S <- init[1]
	D <- init[2]
	p <- params$p
	F <- params$F
	r <- params$r
	beta <- params$discount
	d <- params$d
	m <- params$m
	excess_return <- p - r*F 
	collision_risk_cost <- L(S,D,params)*F 
	congestion_channel <- L_S(S,D,params)*S*F 
	pollution_hazard_channel <- beta*(G_S(S,D,params) + m*(L(S,D,params) + S*L_S(S,D,params)))*L_D(S,D,params)*S*F - (1 - beta)*m*L_D(S,D,params)*S*F
	pollution_persistence_channel <- (1 - d - G_D(S,D,params))*(beta*(p - (L(S,D,params) + L_S(S,D,params)*S)*F) - (1 - beta)*F)
	int_ss <- excess_return - collision_risk_cost - congestion_channel - pollution_hazard_channel - pollution_persistence_channel
	physical_ss_condition <- -d*D + G(S,D,params) + m*L(S,D,params)*S
	return(c(int_ss, physical_ss_condition))
}

# OConstant launch rate steady state finder. Returns an (S,D) pair.
ss_finder_const <- function(init,params,constant_launch_rate) {
	solution <- nleqslv(init,ss_cond_const,params=params, constant_launch_rate=constant_launch_rate)$x
	return(solution)
}

# Open-access steady state finder. Returns an (S,D) pair.
ss_finder <- function(init,params) {
	solution <- nleqslv(init,ss_cond,params=params)$x
	return(solution)
}

# Optimal steady state finder. Returns an (S,D) pair.
opt_ss_finder <- function(init,params) {
	solution <- nleqslv(init,ss_cond_optimal,params=params)$x
	return(solution)
}

# Finds origin point before physical dynamics step. Intended use case: feed it a point on the launch ray leading to the steady state, and it'll undo the physics to show you where the point came from
physics_inverter_objective <- function(pt,S,D,params) {
	S_1 <- pt[1]
	D_1 <- pt[2]
	S_piece <- S_(0,S_1,D_1,params) - S
	D_piece <- D_(0,S_1,D_1,params) - D
	objective <- S_piece^2 + D_piece^2
	return(objective)
}

physics_inverter <- function(pt,params) {
	result <- optim(c(0,0), physics_inverter_objective, S=pt[1], D=pt[2], params=params, method="L-BFGS-B", lower=0, upper=4)$par
	return(result)
}

# Comparative statics inequalities: dX/dS
dXdS_negative <- function(S,D,params) {
	condition <- L_S(S,D,params)*S < 1 - L(S,D,params)
	return(condition)
}

dXdD_negative <- function(S,D,params) {
	X <- eqm_launch_rate(S,D,params)
	S_next <- S_(X,S,D,params)
	D_next <- D_(X,S,D,params)
	condition <- L_S(S_next,D_next,params)*S < ( L_D(S_next,D_next,params)/L_D(S,D,params) )*(1 - d + G(S,D,params))
	return(condition)
}

fancy_Y_LHS <- function(D,params) {
	sats <- eqm_sat_stock(D,params)
	LHS <- d*D
	return(LHS)
}

fancy_Y_RHS <- function(D,params) {
	fe_eqm <- params$fe_eqm
	sats <- eqm_sat_stock(D,params)
	RHS <- G(sats,D,params) + m*fe_eqm*sats
	return(RHS)
}

# analytical expression for fragment autocatalysis rate
FAR <- function(S,D,params) {
	result <- G_D(S,D,params) - d
	return(result)
}


sats_SS <- function(S,D,params,constant_launch=-1) {
	X <- eqm_launch_rate(S,D,params)
	if(constant_launch>0)	{ X <- constant_launch }
	S_imbalance <- S_(X,S,D,params) - S
	return(S_imbalance)
}

debs_SS <- function(S,D,params,constant_launch=-1) {
	X <- eqm_launch_rate(S,D,params)
	if(constant_launch>0)	{ X <- constant_launch }
	D_imbalance <- D_(X,S,D,params) - D
	return(D_imbalance)
}
# analytical expression for kessler syndrome condition

sat_nullcline <- function(D,params) {
	sat_value <- rep(-1,length=length(D))
	for(i in 1:length(sat_value)) {
		result <- uniroot.all(sats_SS,c(0,1e25),D=D[i],params=params)
		if(length(result)==0) {result <- 0} 
		sat_value[i] <- result
	}
	return(sat_value)
}

deb_nullcline <- function(D,params) {
	deb_value <- rep(-1,length=length(D))
	for(i in 1:length(deb_value)) {
		result <- optim(0.1, debs_SS, lower=0, upper=1000, method=c("L-BFGS-B"), D=D[i], params=params)$par
		if(length(result)==0) {result <- 0} 
		deb_value[i] <- result
	}
	return(deb_value)
}
